Cardiac mechanics Benchmark - Problem 3ΒΆ
This code implements problem 3 in the Cardiac Mechanic Benchmark paper
Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S, Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.
[1]:
import dolfin
try:
from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
from dolfin import DirichletBC, Constant, interpolate, Mesh
import pulse
[2]:
geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["benchmark"])
2021-05-13 08:28:25,025 [72707] INFO pulse.utils:
Load mesh from h5
[3]:
# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["CC"] = 2.0
material_parameters["bf"] = 8.0
material_parameters["bfs"] = 4.0
material_parameters["bt"] = 2.0
[4]:
activation = Constant(0.0)
material = pulse.Guccione(
parameters=material_parameters, active_model="active_stress", activation=activation
)
[5]:
# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
return DirichletBC(
V, Constant((0.0, 0.0, 0.0)), geometry.ffun, geometry.markers["BASE"][0]
)
[6]:
# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])
[7]:
# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))
[8]:
# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
[9]:
# Solve problem
pulse.iterate.iterate(problem, lvp, 15.0, initial_number_of_steps=50)
pulse.iterate.iterate(problem, activation, 60.0, initial_number_of_steps=50)
2021-05-13 08:28:25,420 [72707] INFO pulse.utils: Iterating....
2021-05-13 08:28:25,421 [72707] INFO pulse.utils: Current control: 0.000
2021-05-13 08:28:25,421 [72707] INFO pulse.utils: Target: 15.000
2021-05-13 08:28:47,766 [72707] INFO pulse.utils: Iterating....
2021-05-13 08:28:47,767 [72707] INFO pulse.utils: Current control: 0.000
2021-05-13 08:28:47,767 [72707] INFO pulse.utils: Target: 60.000
[9]:
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 294),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 320),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 344),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 370),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 398),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 426),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 454),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 615),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 645),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 804),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 836)],
[Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 292),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 318),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 342),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 368),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 396),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 424),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 452),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 613),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 643),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 802),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 834)])
[10]:
# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)
[11]:
endo_apex_marker = geometry.markers["ENDOPT"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)
[12]:
print(
("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
endo_apex_pos[0]
)
)
Get longitudinal position of endocardial apex: 11.8550 mm
[13]:
epi_apex_marker = geometry.markers["EPIPT"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)
[14]:
print(
("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
epi_apex_pos[0]
)
)
Get longitudinal position of epicardial apex: 15.4904 mm
[15]:
V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
[16]:
from fenics_plotly import plot, set_renderer
set_renderer('notebook')
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()
[ ]: